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Abstract 

Set-membership estimation is usually formulated in the context of set-valued calculus and no probabilistic calculations are necessary. 
In this paper, we show that set-membership estimation can be equivalently formulated in the probabilistic setting by employing sets of 
probability measures. Inference in set-membership estimation is thus carried out by computing expectations with respect to the updated 
set of probability measures V as in the probabilistic case. In particular, it is shown that inference can be performed by solving a particular 
semi-infinite linear programming problem, which is a special case of the truncated moment problem in which only the zero-th order 
moment is known (i.e., the support). By writing the dual of the above semi-infinite linear programming problem, it is shown that, if the 
nonlinearities in the measurement and process equations are polynomial and if the bounding sets for initial state, process and measurement 
noises are described by polynomial inequalities, then an approximation of this semi-infinite linear programming problem can efficiently be 
obtained by using the theory of sum-of-squares polynomial optimization. We then derive a smart greedy procedure to compute a polytopic 
outer-approximation of the true membership-set, by computing the minimum-volume polytope that outer-bounds the set that includes all 
the means computed with respect to V. 
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1 Introduction 

Inferring the value of the state of a dynamical system at the 
various time instants is a classical problem in control and es¬ 
timation theory. The state is estimated based on noisy signal 
observations and on a state transition model, which in turn 
is affected by two sources of uncertainty (namely, process 
disturbance and uncertainty on the initial state conditions). 
In the literature, there are two main approaches for dealing 
with the uncertainties and noises acting on the system: 

• the stochastic (probabilistic) approach that assumes that 
the noises and the uncertainties are unknown but they can 
be described by known probability distributions. 

• the set-membership approach that assumes that the noises 
and the uncertainties are unknown but bounded in some 
compact sets. 

The probabilistic approach is grounded on Bayesian filter¬ 
ing, whose aim is to update with the measurements and 
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propagate up on time the probability density function (PDF) 
of the state. Inferences are then carried out by computing 
expectations with respect to this PDF, i.e., mean, variance, 
credible regions. It is well known that, for linear discrete¬ 
time dynamical systems corrupted by Gaussian noises, the 
Bayesian filter reduces to the Kalman filter. 

The set-membership approach is instead based on the 
construction of a compact set which is guaranteed to in¬ 
clude the state values of the system that are consistent 
with the measured output and the assumed bounds on 
the noises/disturbances [1,2,3,4,5,6]. This compact set is 
propagated in time and updated recursively with the out¬ 
put observations. In set-membership estimation, computing 
inferences thus means to determine this compact set. Set- 
membership estimation was first proposed in [7,8], where 
an ellipsoidal bounding of the state of linear dynamical 
systems is computed. The application of ellipsoidal sets 
to the state estimation problem has also been studied by 
other authors, for example [9,10], and, independently, in the 
communications and signal processing community, starting 
from the works [11,12,13,14]. In order to improve the esti¬ 
mation accuracy, the use of a convex polytope instead of an 
ellipsoid has been proposed in [15,16]. Unfortunately such 
a polytope may be extremely complex and the correspond- 
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ing polytopic updating algorithms may require an excessive 
amount of calculations and storage (without any approxi¬ 
mations, the number of vertices of the polytope increases 
exponentially in time). For this reason, it has been sug¬ 
gested to outer approximate the true poly tope with a simpler 
polytope, i.e. possessing a limited number of vertices or, 
equivalently, faces [17]. In this respect, a parallelotopic 
approximation of the set-membership set was presented in 
[18,19]. A parallelotope is the generalisation of the paral¬ 
lelogram to R". Minimum-volume bounding parallelotopes 
are then used to estimate the state of a discrete-linear dy¬ 
namical system with polynomial complexity. Zonotopes 
have been proposed to reduce the conservativeness of paral¬ 
lelotopes. Intuitively zonotopes are polytopes with parallel 
faces, for a more precise definition see [20, Ch. 2]. A par¬ 
allelotope is thus a special zonotope. Zonotopes are used in 
[21,22,23] to build a state bounding observer in the context 
of linear discrete systems. 

Zonotopes are also employed to address the problem of 
set-membership estimation for non-linear discrete-time sys¬ 
tems with a bounded description of noise and uncertain¬ 
ties [24]. At each sample time, a guaranteed bound of the 
uncertain state trajectory of the system is calculated us¬ 
ing interval arithmetic applied to the nonlinear functions 
through the mean interval extension theorem. This outer 
bound is represented by a zonotope. Similar approaches for 
set-membership estimation for nonlinear systems are pre¬ 
sented in [25,26,27], where ellipsoids are used instead of 
zonotopes. Recently, randomized methods are used in [28] 
to approximate, with probabilistic guarantees, the uncertain 
state trajectory with polynomial sublevel sets. 

The aim of this paper is to address the problem of the es¬ 
timation of the state of a discrete-time non-linear dynami¬ 
cal system (characterized by polynomial non-linearities) in 
which initial state and noises are unknown but bounded by 
some compact sets (defined by polynomial inequalities). We 
are therefore in the context of set-membership estimation, 
but we will address this problem in a very different way 
from the approaches presented above. We reformulate set- 
membership in the probabilistic setting and solve it using 
the theory of moments and positive polynomials. More pre¬ 
cisely the contributions are the following. 

First, by exploiting recent results on filtering with sets of 
probability measures [29,30], we show that set-membership 
estimation can be equivalently formulated in a probabilis¬ 
tic setting by employing sets of probability measures. In 
particular, we show that the prediction and updating steps 
of set-membership estimation can be obtained by applying 
Chapman-Kolmogorov equation and Bayes’ rule point-wise 
to the elements of this set of probability measures V . This 
unifies the probabilistic approach (Bayes filter) and the set- 
membership approach to state estimation. This result can 
have an enormous impact, because it finally can allow us to 
combine set-membership and classical probabilistic uncer¬ 
tainty in order to obtain hybrid filters, i.e., stochastic (prob¬ 
abilistic) filters that are for instance able to use information 


about the bounding region as well as the probabilistic mo¬ 
ments (mean and variance) of the noises or that are able 
to deal with a Gaussian measurement noise and a bounded, 
with known moments, process noise etc.. Moreover, it can 
allow us to compute credible regions (Bayesian confidence 
intervals) that takes into account of both deterministic and 
probabilistic uncertainty, as well as it allows us to make de¬ 
cisions by choosing the action that minimizes the expecta¬ 
tion of some loss function (this is important, for instance, in 
control design). In the context of this paper a first attempt 
in combining deterministic and probabilistic uncertainty has 
been proposed in [29], while [31] has proposed a joint Zono- 
topic and Gaussian Kalman filter for discrete-time LTV sys¬ 
tems simultaneously subject to bounded disturbances and 
Gaussian noises. The work [32] instead proposes a Bayesian 
approach to set-membership estimation imposing a uniform 
distribution on the membership-set similar to the idea pro¬ 
posed in [33,34]. We will show that this approach is differ¬ 
ent from set-membership estimation, since set-membership 
estimation cannot be interpreted in the Bayesian framework, 
but only in the framework of set of probability measures. 
Second, under this probabilistic interpretation, inferences in 
set-membership estimation are carried out by computing ex¬ 
pectations with respect to the set V as in the probabilistic 
case. In particular, we show that the membership set X (i.e., 
the set that includes the state with guarantee) can be obtained 
by computing the union of the supports of the probability 
measures in V. Moreover, we prove that a minimum volume 
convex outer-approximation of X can simply be obtained 
by computing the set Ai that includes all the means com¬ 
puted with respect to the probabilities in V. The proof is not 
constructive, hence we do not have a convenient description 
of Ai. However we show that we can determine the least 
conservative half-space 7 i that includes M, by solving a 
semi-infinite linear programming problem. This problem is 
a special case of the truncated moment problem [35,36,37] 
in which only the zero-th order moment is known (i.e., the 
support). 

Third, by writing the dual of the above semi-infinite linear 
programming problem, we show that, if the nonlinearities in 
the measurement and process equations are polynomial and 
if the bounding sets for initial state, process and measure¬ 
ment noises are described by polynomial inequalities, then 
a feasible solution of the dual can be obtained by simply 
checking the non-negativity of a polynomial on a compact 
set described by polynomial inequalities. An approximation 
of this semi-infinite linear programming problem can be ob¬ 
tained by reformulating it as semidefinite programming by 
using the theory of sum-of-squares (SOS) polynomial opti¬ 
mization. We prove that the approximate solution is robust, 
in the sense that the computed half-space 7 i is guaranteed 
to include Ai, and so the membership set X. 

Fourth, we provide a procedure to determine the minimum- 
volume poly tope S bounding Ai. This procedure is based 
on a refinement of the algorithm originally proposed in [38] 
to compute an approximation of the minimum-volume poly¬ 
tope containing a given semialgebraic set. In particular, we 
use a Monte Carlo integration approach to compute an ap¬ 
proximation of the volume of a polytope, and a greedy pro- 
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cedure to determine an outer-bounding polytope S as the 
intersection of a pre-specified number of half-spaces TLj, 
where each half-space TLj is added to the description of S 
so to minimize the volume of the poly tope including AT 
This allows us to solve the set-membership estimation prob¬ 
lem for polynomial non-linear systems very efficiently and 
through convex optimization. 

Finally, by means of a numerical example involving the 
Lotka Volterra prey-predator model, we show the effective¬ 
ness of our approach. 

2 Problem Description 

Consider an uncertain non-linear discrete-time dynamical 
system described by the difference equations: 

x(fc) = a d (x.(k - 1), k - 1) + w (k - 1), 
y {k) = c d (x(fc),fc) + v(ft), 

where x(fc) = [x±(k),..., a; n (/c)] T G R” is the state of the 
system at the time k, y(k) G R m is the measured output 
vector, w (k — 1) G R” is the process noise and v(k) G 
R" 1 is the measurement noise. In this paper, we consider 
polynomial non-linearities a<j(x(fc), k) and c £ ;(x(fc), k), i.e., 


and therefore n = 2, m = 1, d = 2 and s{d) = = 6. 

We further assume that the only available information about 
the initial state x(0) and the noises w(fc), v(fc) is: 

x(0) G X 0 , w(fc) G Wfc, v(fc) G V fc , (4) 

where Xq, Wfc, Vfc are compact basic semi-algebraic sets, 
i.e., compact sets described by the polynomial inequalities: 

Wfc = {w(fc) G R" : hf(w(k),k) <0, i = 1,... ,t w } , 

(5) 

where hf (with i = 1,... ,t w , t w G N) are polynomial 
functions in the variable w(k). The sets Xq, Vfc are described 
in a similar manner. 

This paper addresses a set-membership filtering problem, 
which aims at recursively estimating, at each time sample 
k = 1,2,, T 0 , (an outer approximation of) the state un¬ 
certainty set Tfc, defined as the set of all values x(fc) com¬ 
patible with the available information, namely the system 
equations (1), the bounds on the initial state and on the noises 
(4), and the output observations y(l), y(2),..., y(T 0 ). For¬ 
mally, the set-membership filtering problem is defined as 
follows. 


a<j(x(fc - 1 ),k - 1) =Afc_iq d (x(fc - 1)), (2a) 

c d {x{k),k) =C fc q d (x(fc)) ! (2b) 

with 

Qd(x) = 

[1,X1, . . . ,X n ,x\,X!X 2 , ■ ■ ■ ,X n -iX n ,X 2 n , 

(3) 

being the vector of all monomials of degrees less than or 
equal to d, which has dimension s(d) = ( n ~^ d ), and Afc_i G 
R” xs ( d ), Cfc G M mxs ( d ) are known time-variant coefficient 
matrices. The resulting system will be referred in the paper 
as uncertain time-variant polynomial system of degree d. 


Example 1 Let us consider the discrete-time polynomial 
system: 

X\(k) = x\{k — 1) (2 — X\(k — 1)) + w\(k — 1), 
x 2 {k) = x\{k - l)x 2 {k - 1) + 0.5a:2(A: — 1) + w 2 {k — 1), 

The output equation is given by: y{k) = X\(k) + x 2 (k) + 
v(k). We can rewrite this system as in (2a)-(2b): 


qd(x) 

Afc_i 


[1,X 1 ,X2,x\,XiX2,x\] T 

0 2 0 -1 0 O’ 

0 0 0.5 0 1 0 


Problem 1 [Set-membership filtering] 

Given the system equations (1), the observations, the bound¬ 
ing sets for the noises Wfc, Vfc and the initial state uncer¬ 
tainty set Xo, compute recursively the state uncertainty set 
44 defined as: 

X k = { x(fc)GR n : x(fc) —a d (x(fc - 1 ),k — l)GWfc_i, 
y (k) - c d(x(fc),fc) G Vfc, 
x(fc - 1) G Afc_i } 


for each k = 1,2,..., T a . ■ 

Note that, in general, the sets X k might be nonconvex and 
their representation can become more and more complicated 
as the time index k increases. Under the assumption that 
X k is bounded, algorithms for computing simple sets (e.g., 
boxes, parallelotopes, zonotops or ellipsoidal regions) outer- 
bounding the state uncertainty sets X k have been then pro¬ 
posed to reduce this complexity. After formulating the set- 
membership filtering problem in a probabilistic setting, this 
paper presents an algorithm for computing (an approxima¬ 
tion of) the minimum-volume polytope outer-bounding the 
sets X k . 

3 A probabilistic framework for set-membership esti¬ 
mation 


Set-membership estimation is usually formulated in the con¬ 
text of set-valued calculus. We will show in the following 
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paragraph that set-membership estimation can be equiva¬ 
lently formulated in the probabilistic setting by employing 
sets of probability measures. Consider the set-membership 
constraint (the time index is dropped for brevity of 

notation) with X C K”. This constraint can be translated 
in a probabilistic setting by saying that the only probabilis¬ 
tic information on the value x of the variable X is that it 
belongs to the set X, or equivalently, 

Pr(X £ X) = Pr(A’) = 1, 

where Pr is a probability measure on X 1 2 3 More precisely 
Pr is a nonnegative Borel measure on X ^ In other words, 
this means that we only know the support of the probability 
measure of the variable X. 

The support does not uniquely define a probability measure, 
as there are an indefinite number of probability measures 
with support T’PI Hence, x € X is equivalent to the con¬ 
straint that the probability measure of x belongs to the set 
Vx(X), that is the set of all probability measures on the 
variable X with support X. Let us define with P the Cumu¬ 
lative Distribution Function (CDF) of the probability mea¬ 
sure Pr. For instance on R we have that P(x) = Pr(—oo, x] 
(this definition can easily be extended to R n ). Then we can 
easily characterize the set of probability measures Vx(X) 
as follows: 

V x (X) = {P: J x dP(x) = 1}, (6) 

where the integral is a Lebesgue-Stieltjes integral with re¬ 
spect to P. Hence, because of the equivalence between Borel 
probability measures and cumulative distributions, hereafter 
we will use interchangeably Pr and P. 

3.1 Inference on the state 

In state estimation, we are interested in making inferences 
about X or, equivalently, computing expectations of real¬ 
valued functions g of X. Since there are an indefinite number 
of probability measures with support X, we cannot compute 
a single expectation of g. However, we can compute upper 
and lower bounds for the expectation of g with respect to the 
probability measures Pr with support X. For instance, the 


1 To clarify this aspect, consider the experiment of rolling a dice. 
Assume that the probability Pr of the outcomes x of the dice 
is completely unknown, then the only knowledge about the ex¬ 
periment is that x £ X = {1, 2,3, 4, 5, 6}, or, equivalently, that 
Pr({l, 2, 3, 4, 5, 6}) = 1. Therefore, the statement Pr(A) = 1 is 
a model for our (epistemic) uncertainty about the probabilities of 
the dice outcomes. We only know that a: € {1, 2, 3,4, 5,6}. 

2 The sample space is R 71 and we are considering the Borel a- 
algebra. X is assumed to be an element of the a- algebra. 

3 The uniform distribution is one of them, but it is not the only 
one. So by considering only the uniform distribution as in [32], 
we loose the full equivalence with set-membership. 


upper bound for the expectation of g is given by the solution 
of the optimization problem: 

SU P Xv5( x ) dP ( x )> 

p (7) 

s.t. p e P x {X), 

which is a semi-infinite linear program, since it has a finite 
number constraints and an infinite dimensional variable (the 
probability measure Pr). Note that we use “sup” instead 
of “max” to indicate that an optimal solution might not be 
attained. The lower bound of the expectation can be obtained 
by replacing sup with inf. 

Problem (7), i.e., determining an upper bound for the expec¬ 
tation of g with respect to the probability measure Pr given 
the knowledge of its support X, is a special case of the trun¬ 
cated moment problem [35,36,37] in which only the zero-th 
order moment is known (i.e., the support). Hence, we have 
the following result [39], [40, Lemma 3.1]: 

Proposition 1 The optimum of (7) is obtained by an atomic 
measuriP] Pr = <5*, where x = argsup xg ^ p(x). 

Note in fact that, VPr £ 'Pa'(X), with associated CDF P, 

E[g} = [ s(x)dP(x) < [ 3 (x)$i(dx) = g(x), 

Jx Jx 

where g(x), by definition of x, is the supremum of g on 
X. The first integral must be understood as a Lebesgue- 
Stieltjes integral with respect to the cumulative distribution 
of an atomic measure on R n . This means that dP(x) denotes 
the distributional derivative of the cumulative distribution 
of an atomic measure, that are in our case Dirac measures 
6x(cbc) (hence the second integral). From this result, it fol¬ 
lows that the probability measures that gives the lower and 
upper bounds for the expectation of g are atomic (discrete) 
measures. 

In order to formulate the set-membership filtering problem 
in a probabilistic framework it is useful to exploit a result de¬ 
rived by Karr in [39], where it is proven that the set of prob¬ 
ability measures VxiX) which are feasible for the semi¬ 
infinite linear program problem (7) is convex and compact 
with respect to the weak* topology. As a result, Vx (X) can 
be expressed as the convex hull of its extreme points and, 
according to Proposition 1, these extreme points are atomic 
measures on X, i.e.: 

V x (X) = Co{S A : x£ X}, (8) 

where = means equivalent in terms of inferences (expecta¬ 
tions). Summing up what we have obtained so far: 


4 An atomic measure in R 71 is a measure which accepts as an 
argument a subset A of R", and returns <5 X (A) = 1 if x £ A, 
zero otherwise. 
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(1) the set-membership constraint x G X is equivalent to 

(6); 

(2) for the inferences, 74 (X) is equivalent to the convex 
hull of all atomic measures on X , (8). 


by applying the Chapman-Kolmogorov equation point-wise 
to the prob ability me asure sPr(-|x(fc — l))in7 : ’(X(fc)|X(fc— 
1)) and Qr in V(X(k — 1)) we obtain 


Hence, we can derive the prediction and updating step 
for set-membership estimation by applying the Chapman- 
Kolmogorov equation and Bayes’ rule to the set of prob¬ 
ability measures in (8). This means that, by reformulating 
set-membership constraints in a probabilistic way, we can 
reformulate set-membership estimation in the realm of 
stochastic (probabilistic) filtering applied to set of probabil¬ 
ity measures. 


Pr(x(fc)) = [ [ / x(fe) (x')dP(x'|x(fc - 1 ))dQ(x(k - 1)) 

jR n JR n 

/ / 4(fc) (x )^ad(x(fc-l),fc-l)+w(^ )<5 x (<ix(fc 1)) 

JR” JR” 

/ ^ad(x(k—i),k— r)+w (x(A;))(5 x [dx(k 1)) 

JR” 


= <5, 


'a<i(x,fc — l)+w 


(x(fc)) 


( 12 ) 


3.2 Propagating in time and updating set of distributions 

We start by deriving the set-membership filtering prediction 
step by applying the Chapman-Kolmogorov equation. 

Theorem 1 (Prediction) Consider the system equation in 
(1) with w (k — 1) G y\4_ i and assume that the only prob¬ 
abilistic knowledge about X(fc — 1) is the support X k _\. 
Then it follows that the probability measure Pr on the value 
x(k) of the state at time k belongs to the set 

^(X(fe)) = Co{4: xG^}, (9) 

with 

X k = jx(fc) : x(k) = a d (x(k - 1), k - 1) + w (k - 1) 

with x(k — 1) G X k ~i 1 w (k — 1) G W4-i j, 

( 10 ) 

or equivalently: 

X k = jx(fc) : x(fc) - a d (x(k — 1), k — 1) G W k -\ 

with x(k — 1) G (11) 

Proof: Let us consider the time instant k. From the system 
equation in (1), w(fc — 1) G VV4_i and (8), it follows that 

V(X(k)\x(k-l)) 

= 1 ‘ {^ad(x(fc — l),k— l)-{-w • ^ C Wfc—l} , 

this is the conditional set of probability measures for the vari¬ 
able X(fc) given the value x(k — 1) of the variable X(fc — 1) 
Hence, since X(fc — 1) G X k _i and so the set of probability 
measures for the variable X(fc — 1) is 

Vx^X.ik- 1)) =Co{5 r : x G **_!}, 


where I^{k) (x r ) denotes the indicator functiorFl and with 
x G X k _i and w G W4 _i and where we have ex¬ 
ploited the fact that f Rn 4(fe) (x , )4 d (x(fc-l),fe-l)+w(rfx , ) = 

4(fc)(aj(x(fc-l),fc-l) + w) = <5 ad (x(fe-l),fe-l)+w(x(fc)). 

From (8), (12) and the definition of X k , the theorem fol¬ 
lows. ■ 

Theorem 1 shows that, by applying the Chapman- 
Kolmogorov equation point-wise to the probability mea¬ 
sures in P(X(k)\x(k — 1)) and 7 ? Xk _ 1 (X(k — 1)), we can 
obtain a set of probability measures Vg (X(fc)), which is 
completely defined by its support and whose support coin¬ 
cides with the one obtained in set-membership estimation 
after the prediction step. 

We now derive a similar result for the updating step. 

Theorem 2 (Updating) Consider the measurement equa¬ 
tion in (1) with v(fc) G 14 and assume that the only proba¬ 
bilistic knowledge about x(fc) is described by (9)-(10). Then 
it follows that the updated probability measures P on the 
value x(k) of the state at time k belongs to the set: 

V Xk (X{k)) = Co{8^-. xG^}, (13) 

where 

x k = x k cy k , (14) 

with 

34 = {x(fc) : y(fc) - c d (x(k),k) G 14}. (15) 


Proof: Observe that, at each time /;:, 

P(Y(fc)|x(fc)) = Co {8 Cd(x(fc) , fc)+ v : v G 14} . 

Then, the updating step consists of applying Bayes’rule 
to the probability measures V(Y(k)\x(k)) and to Q in 


5 fx(fc)(x / ) = 1 when x(fc) = x' and zero otherwise. 
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by X k , i.e.. 


(17) 


dP(x(k)\y(k)) 


fu™ Iy(k)(y')dP(y'\x(k))dQ(x(k)) 
fun f R m Iy(k){y')dP(.y'\x(k))dQ(x(k)) 
Pr(y(fc)|x(fc))dQ(x(fc)) 
fun Pr(y(k)\x(k))dQ(x(k)) ’ 


where we have exploited the fact that 


Pr(y(fc)|x(fc)) = f Iy(k) {y')dP(y l \x(k)). 

J R m 


Note that the probability of a point on R" can be nonzero 
since Pr is an atomic measure. In order to apply Bayes’ rule 
we need to ensure that the denominator is strictly greater 
than zero: 


fun Pr(y(k)\x(k))dQ(x(k)) 

= fu n dc d (x(k),k)+v(y(k))dx(dx(k )) > 0 . 



dP(x(k)) = 1, 


where X k is given by (14), or equivalently by (6). In other 
words, the support of the probability measure Pr of the value 
of the state x(k) given the output observation y(k) and the 
system equations (1) is nothing but X k . This is in accordance 
with the set-membership formulation, which claims that 
x(/c) belongs to state uncertainty set X k defined in (6). Then 
we can solve set-membership filtering by applying recur¬ 
sively Theorems 1 and 2, as described in Algorithm 1. The 


Algorithm 1: prediction and updating 

Al.l Initialize "P^ o (X(0)) = Co{S R ■ x G X 0 }. 

A1.2 For k = 1, ... ,T 0 : 

Al.2.1 P^ k (X(k)) = Co {<5* : x e A’fej with X k de¬ 
fined in (10); 

Al.2.2 V Xk (X{k)) = Co {6a : x G X k } with X k de¬ 
fined in (14). 


Hence, the above inequality holds if and only if x and v are 
chosen, at time k, such that: 

c d (x, k) + v = y(k). (16) 

Bayes’ rule is only defined for those probability measures 
for which the denominator is strictly positive, that implies 
that the above equality must be satisfied 01] The equality (16) 
can be satisfied only if x G 34 which, together with the 
constraint x G X k , implies that 

x G x k n34- 

Under the constraint (16), it follows that 6 Cd ^.k)+v(y{k)) = 
1 and, thus, the denominator is equal to one. Hence, we have 
that 

dP(x(fc)|y(fc)) = f Rm I y ( k )(y')dP(y'\x(k))dQ(x(k)) 

= fu m dy(k)(y )6 Cd (x(k),k)+v(y(k))Sx(dx(k)) 
6 Cd (x(k),k)+v(y{k))Sa(dx(k)) 
dc d (x,k)+v(y {k'))Sa(dx(k')') 

= 6x(dx(k)) 

with x G X k H y k - Hence, the updated probability mea¬ 
sure Pr(-|y(fc)) on the values of the state at time k is 
Pr(-|y(fc)) = da, which proves the theorem. ■ 

From Theorem 2, the support of the updated probability 
measure Pr on the value x(fc) of the state at time k is given 


This way of updating set of probability measures has been 
proposed by Walley [41, Appendix J] under the name of regular 
extension. 


steps Al.2.1 and Al.2.1 are the prediction and the updating 
steps, respectively. Note that the set of probability measures 
Vx k (X(fc)) (or Vx k (X(fc))) i s computed by taking into ac¬ 
count all the observations y k = {y(l), y(2),..., y(fc)} (re¬ 
spectively y fc_1 ). Hence, it should be more correctly denoted 
as Vx k (X(k)\y k ) (respectively Vx k (X(fc)|y fc_1 )). We have 
omitted this notation for brevity. 

Remark 1 Under the assumptions (2a), (2b) and (5), the set 
Xk is a semialgebraic set in R™, described by the intersec¬ 
tions of the semialgebraic sets X k (Eq. (11)) and 34 (Eq. 
(15)). Formally, X k is the projection in the space ofx(k) of 
the set 

X k = {x G M 2 " : h s (x(k)) <0, s = 1,..., m} , (18) 

where x(k ) is the augmented state vector x(k) = 
[x T (fc) x T {k — 1)] T an d h s (x(k)) (with s = 
are the polynomial functions in x{k) and x(k — 1) (or equiv¬ 
alently in x(k)) defining X k and yk- In the rest of the paper, 
we will use the following notation to describe the set X k : 

X k = (x(fc) G M" : h s {x(k)) <0, s = 1,..., m} . (19) 

Remark 2 The reformulation of set-membership in the 
probabilistic framework is important for two main reasons. 
First, it allows us to reinterpret the operations performed 
in set-membership estimation and justifies them in terms 
of a probabilistic framework. We have just seen the rein¬ 
terpretation of prediction and updating in terms of the 
Chapman-Kolmogorov equation and Bayes’ rule. We will 
further investigate this interpretation in the next sections. 
In particular, in Section 4, we will show that the convex 
membership set computed in set-membership estimation can 
also be interpreted as the set of posterior means calculated 
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with respect to the posterior set of probability measures 
Vx k (X-(k)) (in the Bayesian setting, we know that the 
posterior mean is the optimal estimate with respect to a 
quadratic loss function - a similar result holds for the set 
of posterior means [42, Sec.5]). This result can also now 
be applied to set-membership estimation because, after this 
probabilistic interpretation, we are now able to compute 
expectations. Moreover, in Section 5 we will also highlight 
the connection between set-membership estimation and the 
theory of moments (through duality). 

Second, we are now potentially able to combine set- 
membership and classical probabilistic uncertainty in order 
to obtain hybrid filters, i.e., stochastic (probabilistic) fil¬ 
ters that are for instance able to use information about 
the bounding region as well as the probabilistic moments 
(mean and variance) of the noises or that are able to deal 
with a Gaussian measurement noise and a bounded, with 
known moments, process noise etc.. A first attempt in this 
direction is described in [29] for scalar systems. We plan to 
further investigate this direction in future work by using the 
theory of SOS polynomial optimization, that we also use in 
the next sections. 

4 Computing the support as an inference on the set of 
probability measures 

In the probabilistic formulation of filtering, all the available 
information at time k is encoded in the posterior probability 
distribution of the state x(k) given all the observations y fe ). 
In the set-membership setting, this information is encoded in 
the updated set of probability measures Vx k (X(fc)). Infer¬ 
ences can then be expressed in terms of expectations com¬ 
puted with respect to this set. The set-membership estima¬ 
tion problem can, for instance, be reformulated as follows: 


Theorem 3 Assume that A')- is compact and that Q | C M n 
is a convex set defined as follows: 

fli = arg inf f dx(k ) 

£~2CR n ,f2 conv. q 

s.t. (21) 

fdP(x(k)) = l, VP G Px k (X(k)) 
n 

Then, it results that fli = A i, with 

M = \ J x{k)dP(x{k )): P€Px k (X(fc)) l . (22) 


Proof: From (21) it follows that 0 | is the minimum volume 
convex set that includes X k . Thus, if A'/, : is convex, then 
Qi = Xk- Hence, from (8), the equality 

J x(fc)4 (k )(dx(fc)) = x(fc), 

x k 

and (22), it immediately follows that AA = fl 1 . Conversely 
assume that Xk is not convex, then 1 1 \ D A). . Since O j 
is the minimum volume convex set that includes Xk, then 
Qi must be equal to the convex-hull of Xk- This means 
that for each x £ fii, there exist zi,Z 2 £ Xk such that 
wz\ + (1 — w) Z 2 = x for some w £ [0,1] (by definition of 
convex hull). Then, consider the probability measure 

w8 Zl + (1 — w)d Z2 . (23) 


n* = arg min J dx(k) 
ncR- n 

s.t. (20) 

fdP(x(k)) = 1 , VPeVxMk))- 

n 

The solution of (20) is the minimum-volume set f l C R n , 
such that Pr(x(fc) £ Q) = 1 for all probability measures Pr 
in Px k (Mk)) (i.e., with support TVlPlThus. f l* coincides 
with Xk- Since X k may be not convex, the problem (20) is 
in general difficult to solve. However, the problem can be 
simplified by restricting Q to be convex, thus computing a 
convex outer-approximation of Xk- 

The following theorem shows that computing the minimum- 
volume convex set Q such that P(x(fc) £ fl) = 1 is 
equivalent to obtain the set that includes all the possible 
means computed with respect to the probability measure in 

Vx k (Mk)). 


' It is thus the union of all the supports of the probability measures 

in Px k (X(fc)). 


Because of (8), it holds: 

w6 Zl + (1 - w)S Z2 £ Vx h (Mk)), (24) 
and 

J x(k) (■wS zl (dx(k )) + (1 - w) 5 Z2 (dx(k))) = x. (25) 
x k 

Thus, x belongs to M, and vice versa. ■ 

Theorem 3 has the following fundamental implications: 

• a convex outer-bounding of the set of all the possible 
means computed with respect to the probability measures 
in TAtfc(X(fc)) (i.e., the set M) is also a convex outer- 
bounding of the support Xk of the set of probability mea¬ 
sures V Xk (X.(k)). 

• the tightest convex outer-bounding of the support Xk of 
the set of probability distributions 7^ (X(fc)) is the set 
of the means computed with respect to the probability 
measure in 7 ? ^ i .(X(fc)). 
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We can thus use A4 as an outer-approximation of A).. Al¬ 
gorithm 1 is therefore modified to include the following ad¬ 
ditional steps. 


Refinement of Algorithm 1: outer-approximation step 

Al.1.3 Outer-approximate A'’;,, with M defined in (22). 
Al.1.4 Redefine Vx k (X(k)) = Co{Sx : x G M}. 


M (or equivalently, including X), is obtained for v = v*, 
with 

v* = max / ut T xdP(x) 

s.t. f dP(x) = 1. (28) 

Proof: Let p = f xdP(x) be a point belonging to A4. Let 
us first prove that if v > v*, then M. C 77. First, note that: 


Unfortunately, Theorem 3 does not provide a constructive 
way to find the set M. However, by restricting the outer- 
approximation of the support Xk to have a simple form (e.g., 
a polytope). Theorem 3 can be still exploited to determine an 
outer-bounding set of Xk- The following theorem provides 
results to compute an outer-bounding box of Afc. 

Theorem 4 (Box approximation) The minimum volume 
box that includes Xk can be found by solving the following 
family of optimization problems 


ui T p < v* =supu> T / xdP(x) 

p J 

s.t. f dP(x) = 1 
Jx 

Therefore, for v > is*, u> p < is* < is, which means that 
p = f xdP(x) also belongs to 77 for all p G M. Thus, TL 
contains AL By choosing is = v*, we obtain the tightest 
half-space defined by the normal vector uj that includes AL 


x*(k) = opt f Xi(k)dP(x(k )) 

. P (26) 
s.t. f dP(x.(k)) = 1. 

for i = 1,... ,n, where by selecting opt to be min or max 
we obtain the half-spaces J Xi(k)dP(x(k )) > x*(k) and, 
respectively, f Xi(k)dP(x(k)) < x*(k) which define the 
box. 

The proof of Theorem 4 is provided together with the proof 
of Theorem 5. Based on Theorem 4, by computing the lower 
and upper means of the components x\ (fc),... , x n (k) of the 
vector x(k), the tightest box that outer-approximates Xk is 
obtained. In the following we will discuss how to efficiently 
solve optimization problems similar to (26) and how to find 
an outer-approximation of Xk that is less conservative than 
a box. For simplicity of notation, in the rest of the paper, the 
dependence of the state x(fc) and of the set A), on the time 
index k will be dropped, and only used when necessary. 


5 Exploiting duality 

In this section we discuss how to efficiently solve optimiza¬ 
tion problems similar to (26). In particular, we slightly mod¬ 
ify (26) in order to be able to determine the more general 
half-space 

U = {p G R n : w T p < u} , (27) 

where ix> G R”, ^Gl and p = f xdP(x )0 

Theorem 5 Let us fix the normal vector u> defining the half¬ 
space Tl in (27). Then, the tightest half-space IT including 


It can be observed that (28) reduces to (26) when uj = e* 
for i = 1,..., n, where e, is an element of the natural 
basis of R". Note that, in Problem (28): (i) the optimization 
variables are the amount of non-negative mass assigned to 
each point x in X (i.e., the measure Pr(x)); (ii) the objective 
function and the constraint are linear in the optimization 
variables. Therefore, (28) is a semi-infinite linear program 
(i.e., infinite number of decision variables but finite number 
of constraints). By exploiting duality of semi-infinite linear 
program (see for instance [43]), we can write the dual of 
(28), which is defined as: 

is* = inf is 

V 

s.t. is > tu T x, Vx G X , (29) 

which is also a semi-infinite linear program (i.e., finite num¬ 
ber of decision variables (is) but infinite number of con¬ 
straints). A solution is is feasible for Problem (29) provided 
that: 

is — tu T x >0, Vx G X. 

Hence, checking the feasibility of is is equivalent to check 
the non-negativity of the polynomial is —u> T x in the set A". 

Remark 3 The probabilistic formulation of the set- 
membership estimation described so far is general enough, 
and it is valid also when the dynamical system in (1) is 
not a polynomial system and when the uncertainty sets 
Ab,Wfc,Vfe in (4) are not semialgebraic, but just com¬ 
pact sets. The assumptions of polynomiality are used in 
the following to efficiently solve the semi-infinite linear 
programming problem (29) through convex optimization. 

5.1 Sum-of-squares polynomials 

A sufficient condition for a polynomial to be non-negative 
over a semialgebraic set is that it can be written in terms of 
sum-of-squares (SOS) polynomials (see, e.g., [44]). 


The half-space IT lies on the space of the means. 
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Definition 1 A polynomial <r(x), with x g R 2 ", of degree 
2d is a sum-of-squares polynomial, denoted by er(x) g £[x], 
if and only if it can be written as: 

cr(x) = q d (x) T Qq d (x), (30) 

where Q is a real symmetric positive semidefinite matrix 
of dimension ( 2ra ,^ d )- The vector of monomials q d (x) is 
defined as in (3). The set of SOS polynomials of degree less 
then or equal to 2d is denoted as £ 2d [x]. 

Then, for a given integer d > 1, a sufficient condition for 
is — u; T x to be non-negative in X is (see for instance [37, 
Ch. 4]): 

m 

v — ut T x = ct o(x) — V a s (x.)h s (x) V x g R 2n 

s=i (31) 

CT 0 (x),o-i(x),... ,er m (x) g £ 2d [x], 

where h s (x.) (with s = 1,..., to) are the polynomial non¬ 
positive inequality constraints defining the semialgebraic set 
X. In order to avoid confusion, we would like to stress 
that also v — u; T x is a polynomial in the variable x. In 
fact, we remind that the augmented state x(k) is defined as: 

x(fc) = [x T (fc) X T (fc — 1)] T - 

The following (more conservative) optimization problem can 
be then solved instead of (29): 

v** = inf i/ 

V,G S 

m 

V — u; T X = CTo(x) — <T s (x)/l s (x), V X g l 2 " 

6 — 1 

CT 0 (x),CTi(x), . . . ,CT m (x) g £ 2d [x], 

(32) 

Note that, by rewriting the SOS polynomials <r s (x) (with 
s = 0,..., to) as in (30), Problem (32) can be also rewritten 
as: 

v** = inf v 

v - w T x =q d (x) T Q 0 q d (x)+ 

m 

~y^qd(x) T Q s q d (x)/t s (x), Vxgl 2 " 

6=1 

Q s A 0, s = 0,..., m. 

(33) 


Some remarks: 

(1) Problem (33) is a semidefinite programming (SDP) 
problem [44,45], thus convex. In fact, checking if the 
polynomial is — u> T x is equal to q d (x) T Qoq d (x) — 
£6 m =iq d (*) T Q6q d (x)Mx) for all x g R 2 " leads 


to linear equalities in is and in the matrix coeffi¬ 
cients Q s (with s = 1, ...,m). Besides, enforcing 
cto(x), oi(x),..., cr m (x) to be sum of square poly¬ 
nomials leads to linear matrix inequality (LMI) con¬ 
straints in the coefficients of cto(x), <ti(x), ..., a m {x) 
(i.e., Q s >z 0). 

(2) For v = is**, the robust constraint is** — w T x > 
0 Vx g X appearing in Problem (29) is guaranteed to 
be satisfied. As matter of fact, for all x g X, h s (x) < 0 
(with s = 1,..., m) by definition of X. Furthermore, 
the SOS polynomials cr s (x) = q d (x) T Q s q d (x) (with 
s = 0,... ,m) are always nonnegative over R 2ra as 
Q s A 0. Thus, both the left and the right side of the 
equation in Problem (33) are nonnegative for all x g X. 

(3) Since the equality constraint in (33) gives only a suffi¬ 
cient condition for the non-negativity of is — ui' x on 
X, it follows that is* < is**. Therefore, conservative¬ 
ness is introduced in solving (33) instead of (29), as 
highlighted in Corollary 1. 

(4) However, according to the Putinar’s Positivstellensatz 
(see, e.g., [46] and [47, Ch. 3]), a polynomial which is 
nonnegative over a compact semialgebraic set X can 
exactly always be written as a combination of SOS 
polynomials, provided that the degree of the SOS poly¬ 
nomials ct o (x),..., a m (x) is large enough. In other 
words, we can make is** close to is* by increasing the 
degree of the SOS. However, in practice it often hap¬ 
pens that the relaxed solution is** and the optimal one 
is* coincide with each other for small values of the SOS 
degree 2d. 

Corollary 1 The set M. is guaranteed to belong to the half¬ 
space TT : u; T x < is**, i.e. 

MCJi. (34) 


Proof: The proof straightforwardly follows from Theorem 5 
and is* < is**. U 

Example 2 Let us consider the discrete-time polynomial 
system described by the difference equations: 

x\{k)=x\(k—l)x2(k—l\xi{k—T) + x 2 (k — l))+wi(k-l), 

x 2 {k)=xi(k—l)x2{k—l\2x\(k—T) + x-fk-l))+w 2 (k— 1) 

(35) 

The output equation is given by: y(fc) = X\(k) + x 2 {k) + 
v(fc). The following conditions are assumed: (i) the ini¬ 
tial state x(0) belongs to Xq = {x(0) : Jlx(0)|| 2 < 0.2}, 
the process noise w (fc) = [iui(fc) w 2 (k)]' is bounded by 
l|w(fe)|| 2 < 0.4, and the measurement noise by ||v(fc)|| 00 < 
0.5. The observed output y (k) at time k ~ 1 is y(fc) = 0. 
We are interested in computing an half-space PI : os' p < is 
containing the state uncertainty set X^ (or equivalently JVl) 
at time k = 1. The normal vector u> characterizing PI is 
fixed and it is equal to uj = [— 1 — 0.5] T . In order to com¬ 
pute the constant parameter is defining Pi, the SDP Problem 
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Fig. 1. True state uncertainty set X\ (dark grey region) and half-s¬ 
pace PI : —pi — 0.5p2 < 0.45 (light gray region). 

(33) with x(l) = [x T (l) x T (0)] T and 

h\ (x(l)) :xi(0) 2 + x 2 (0) 2 - 0.2 2 < 0, (36) 

L 2 (x(l)):(xi(l)-xi(0)x 2 (0)(xi(0) + x 2 (0))) 2 + 

'-V-' 

»i 2 (0) 

(®i(l)-®i(0)a: 2 (0X2a;i(0) + x 2 (0))) 2 - 0.4 2 <0, 
^---' 

«>i(o) 

hs (x(l)) :y( 1) - aci(l) - x 2 (l) - 0.5 < 0, 

v(l) 

ft 4 (x(l)): - ( y(l) - xi(l) - x 2 (l) ^ - 0.5 < 0, 
v(l) 

(37) 

is solved for a SOS degree 2d = 4. The SOStools [48] has 
been used to easily handle the SOS polynomials appearing 
in (33). The CPU time taken by the solver SeDuMi [49] 
to compute a solution of the SDP Problem (33) on a 2.40- 
GHz Intel Pentium TV with 3 GB of RAM is 2.1 seconds. 
The computed half-space TL is plotted in Fig. 1, along with 
the true state uncertainty set X\. According to Theorem 5 
and Corollary 1, X\ is included in the half-space TL. Note 
also that, although the original robust optimization problem 
(29) has been replaced with the SDP problem (33), the com¬ 
puted parameter it** defining FI is such that the hyperplane 
cu T x = v** is “almost” tangent to the set X\. Thus, only a 
small level of conservativeness is introduced in using SOS. 

6 Computation of the minimum-volume polytope con¬ 
taining M. 

In the previous section, given the normal vector ut defining 
the half-space PL in (27), we have shown how to compute, 
through convex optimization, the constant parameter v such 
that M C PL. 


Now consider the following family of half-spaces: 

PLj = {p € R" : utjp < vj} , 

for j = 1, J with J > n + 1. Our goal is to choose the 
normal vectors utj, along with the constant parameters v 3 , 
defining the half-spaces PLj such that 

(1) MCS = n / =1 Hf, 

(2) the polytope S has minimum volume. 

In other words, now also the normal vectors utj for j = 
1,..., J have to be optimized. Then, we can formulate the 
problem we aim to solve as: 

inf / g?x s.t. M C S, (38) 

J s 

where S in (38) is constrained to be a polytope. There are 
two main aspects making (38) a challenging problem, i.e., 

(1) the minimum-volume poly tope outer-approximating a 
generic compact set in M" might not exist. For instance, 
if A4 is an ellipsoid, its convex hull is described by an 
infinite number of half-spaces, namely all the support¬ 
ing hyperplanes at every boundary point of A4. 

(2) the problem of computing the exact volume J s r/x of 
a polytope S in R" is fiP -hard (see, e.g. [50,51]. The 
interested reader is also referred to [52] for details on 
^P-hard problems). Although several algorithms have 
been proposed in the literature to compute the volume 
of a polytope S through triangulation [53,54,55,56], 
Gram’s relation [57], Laplace transform [58] or ran¬ 
domized methods [59,60,61], all the approaches men¬ 
tioned above require an exact description of the poly¬ 
tope S in terms of its half-space or vertex representa¬ 
tion. However, in our case, the parameters utj , Vj defin¬ 
ing the half-spaces PLj are unknown, as determining 
utj , itj is part of the problem itself. 

In the following paragraph we present a greedy algorithm to 
evaluate an approximation of the minimum-volume poly tope 
outer-approximating the set A4. 

6.1 Approximation of the objective function 

As already pointed out in the previous paragraph, one of the 
main problems in solving (38) is that an analytical expression 
for the computation of the volume of a polytope S in JR" is 
not available and the polytope S is unknown, as computing 
S is part of the problem itself. In order to overcome such 
a problem, a Monte Carlo integration approach [62] is used 
here to approximate the volume of S. Specifically, given an 
outer-bounding box B of the set A4 (which can be computed 
as discussed in Theorem 4) and a sequence of N random 
points independent and uniformly distributed in B, 
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the integral f s dx can be approximated as: 


N 

/ dx ps Vol(B) — (39) 

-' 5 V *=1 


where Vol(B) is the volume of the box B and I(pf) is the 
indicator function of the (unknown) polytope S defined as 


Algorithm 2: Polytopic outer approximation S of M 
[input] List C = {pi}fLi of N random points uniformly 
distributed in the box B. 

A2.1 Set j = 1. 

A2.2 Compute the half-space PLj, defined as PLj : ujJ p — 
Vj < 0 (with u>j ^ 0), that contains the minimum number 
of points in the list £ and such that A4 is included in PL j, 
i.e., 


J {S}(Pi ) 


1 if Pi £ S 
0 otherwise 


(40) 


Remark 4 It is worth remarking that: 


E 


1 

V i =1 


Vol{S), 


where the expectation is taken with respect to the random 
variable p^ Furthermore, because of the strong law of large 
numbers, 


N 


lim Vol(B)—y^I {s} {pi) = Vol(S) w.p. 1, (41) 

N—foo N z ' L J 


N 


i= 1 


N 

u i > v *o =ar 8 ““ Z hn^iPi) 

OJ j tlK 

^■eR 1=1 

s.t. (43) 

ujj 0 

M C PLj 

Pi € £, i = l,...,N 

A2.3 Collect all the points pi £ £ belonging to the half¬ 
space PLj (computed through (43)) in a list £j. Let Nj be 
the number of elements of £ t . 

A2.4 If Nj < N, then £ <— £j, N t— Nj, j •<— j + 1 and 
go to step A2.2. Otherwise, set J = j — 1 and go to step 
A2.5. 

A2.5 Define the polytope S as S = B f~l Dj^i Hj- 

[output ] Polytope S. 


where w.p. 1 is for with probability 1. For finite samples N, 
the level of accuracy of the approximation in (39) depends 
on the shape of the set S as well as on the volume of the 
outer box B. The reader is referred to as [62] for details on 
Monte Carlo integration methods. 

On the basis of (39), the volume minimization of problem 
(38) can be then approximated as 

N 

min^/{ S }(pi) s.t. M C S (42) 

i=1 


In the following subsection, we describe a greedy procedure 
aiming at computing an approximation of the minimization 
problem (42). 

6.2 A greedy approach for solving (42) 

The key steps of the approach proposed in this section to 
compute a polytopic outer-approximation S of the set M 
are summarized in Algorithm 2. 

Algorithm 2 generates a sequence of half-spaces PLi,..., PLj 
as follows. First, the half-space PL\ that minimizes an 
approximation of the volume of the poly tope B FI PL\ is 
computed. The approximation is due to the fact that the 
volume of B Cl PL\, given by the integral dx, is ap¬ 

proximated (up to the constant Vo1 ^ ) by J2iLi I {Hi} {‘Pi) 


(corresponding to the objective function of problem (43)). 
Then, the new half-space PL -2 that minimizes an approxima¬ 
tion of the volume of the polytope BnPLi ITH2 is generated. 
In order to approximate the volume of B Cl PL 1 (T PL 2 , all 
the points pi of the list £ = that do not belong 

to the polytope B IT PL\ are discarded, and all and only 
the points belonging to B IT PL 1 are collected in a new list 
£\ = {pi}^!fi (step A2.3). The volume of B H PL\ D PL 2 

is then approximated by I{H 2 }(Pi)’ Pi £ £ 1 . 

The procedure is repeated until iVj + i = Nj (step A2.4), 
which means that the number of samples pi belonging to 
the polytope B f~l PLi IT ... fl PLj+ 1 is equal to the number 
of samples p, belonging to the polytope B fl PLi fl... fl PLj. 
Note that, because of the constraint A4 C PLj appearing in 
optimization problem (50), the half-spaces PL \,..., PLj are 
guaranteed to contain the set A i, and thus S = Bflfj ; J =1 PLj 
is an outer approximation of AL Finally, we would like 
to remark that, in case we are interested also in bounding 
the maximum number of half-spaces defining the polytopic 
outer approximation S, Algorithm 2 can be stopped after 
an a-priori specified number of iterations. 

Example 3 Let us consider again Example 2. The first steps 
of Algorithm 2 are visualized in Fig. 2. An outer-bounding 
box B of the true state uncertainty set (dark gray region) is 
first computed (Fig. (a)). A set of 80 random points (black 
dots) uniformly distributed in B is generated (Fig. (b)). The 
half-space PL\ containing the true state uncertainty set and 
the minimum number of points is computed. The points which 
do not belong to PL\ are discarded (gray dots in Fig. (c)). A 
new half-space PL 2 containing the true state uncertainty set 
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and the minimum number of black dots is computed (Fig. 
(d)). Again, the points that do not belong to PL\ fl PL 2 are 
discarded (gray dots in Fig. (d)). The procedure terminates 
when no more black points can be discarded. 
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(a) 
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(b) 



(c) (d) 


Fig. 2. First steps of Algorithm 2. 



Fig. 3. Indicator function I-h- ( pi ) (black solid line) and approxi¬ 
mate function R{n }(pi) (gray thin line). When uijpi — Uj > 0, 
I{Hj}(Pi) an d are overlapped and they are equal to 0. 

tions I{-u.}(pi) with the convex functions R^ j y(pi), i.e., 

N 

> ** = al 'g min 51 R {Hi } ( Pi ) 

Wjtl 
i/jgi 1=1 

s.t. (46) 

ujj ^ 0 

M C Uj 

Pi e A * = 1 ,...,7V. 

Theorem 6 If(i) there exists at least one point pi in the list 
C such that uj* pi — u* <0 and (ii) uj* , u* is the optimal 

solution of problem (46), then the hyperplane uj* p—v* =0 
is a supporting hyperplane for the set A4. 


Technical details of step A2.2, which is the core of Algo¬ 
rithm 2, are provided in the following sections. 


6.3 Approximation of the indicator functions 


Note that the objective function of problem (43) is noncon- 
tinuous and nonconvex since it is the sum of the indicator 
functions I{ H .}(pi) defined as 


hu^Pi) 


1 if ujjpi - Uj < 0, 
0 if ujpi - Uj > 0. 


(44) 


We then transform it in a convex objective function. Each 
indicator function Isj ij ^(pi) is here approximated by the 
convex function Rp H . j(xi) defined as 




-ujJ Pi + Uj if ljJ Pi - Uj < 0 , 
0 if ajjpi - Uj > 0. 


(45) 


A plot of the functions I^ }(pi) and Rpu } (Pi) is given in 
Fig. 3. 


Problem (43) is thus relaxed by replacing the indicator func- 


Proof: Theorem 6 is proved by contradiction. Let H* be the 

half-space defined as PL* : uj* p — P*< 0. Let us suppose 
that uj* ,Uj is a feasible solution of problem (46) such that 

uj* p — u* = 0 is not a supporting hyperplane for A4, that 

is, for some £ > 0 ,PLj : uj* p—u*+e< 0 for all x £ A4. 
Let us define Uj as Uj = u* — e. Note that {uj*,i>j} is 
still a feasible solution of problem (46) and PLj C PL*. Let 
V* = Y^iLi be the value of the cost function of 

Problem (46) obtained for uj = uj* and u = u*. 
is then given by 




-ujfpi + u* if ujfpi - u* < 0 
0 if uj* pi - u* > 0 


(47) 


Similarly, let V = Y^=iR{H }ijPi) be the value of the 
cost function of Problem (46) obtained when uj = uj* and 
u = Uj. The term R^p ^pi) is the given by 


R 


{Hi} 


[0 if uj* pi — Uj > 0 


Since PLj C PL*, then when R^*y(pi) = 0, also R^.y(Pi) 
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is equal to zero. On the other hand, when = 

—u)* pi + v* > 0, then R^-p \('Pi) can be equal either to 

zero or to —Co* pt+bj = -lj* Pi + v*-e <-lj* Pi+v*. 
On the basis of the above considerations, it follows: 

f R { * ;} (Pi) = ^ Pi - v*j > 0 (49) 

\ > R {H 3 }(ft) if Q f Pi ~ D j < 0 

Since by hypothesis (i) there exists at least one point pi in 
the list C such that Co* Pi~v* < 0, it follows that V* > V. 
Therefore, Co* , v* is not the optimal solution of problem 
(46). This contradicts hypothesis (ii). ■ 


following SDP problems: 

N 

To *, V* = arg jmm n ^ R {Uj } (ft) 

lAjER 1=1 

Qs 

s.t. 

LO jt l = 1 

Vj - w,,x = q d (x) T Qoq d (x)+ 

m 

-^ q d(x) T Q s q d (x)/i s (x), Vx€l 2 " 

.9=1 

Q S h 0, s = 0,...,m. 

Pi G C, i = l,...,N 


(51a) 


Theorem 6 has the following interpretation. Among all 
the half-spaces defined by the normal vector to* and con¬ 
taining the set A4 , the optimization problem (46) provides 

T 

the half-space PL* : Co ■ p — u* < 0 which minimizes the 
volume of the poly tope B fl 'H\ D ... IT 'H* , even if the 
integral f Bnn * n dx is approximated (up to a constant) 

by e;= ■iI{W}(Pi) anc * the indicator functions I{u*}(Pi) 
are replaced by the convex functions 7?{-H*}(ft). 


6.4 Handling the constraint A4 C R :j 


The constraints A4 C PLj can be handled through the SOS- 
based approach already discussed in Section 5.1. Specifi¬ 
cally, by introducing a SOS relaxation, Problem (46) is re¬ 
placed by: 


N 

to* , v* = arg min R {H] } (p« ) 

wj-eR 1=1 

Qs 

s.t. 

LOj ^ 0 

T (50) 

Vj - WjX = qd(x) Qoq d (x)+ 

m 

-^q d (i) T Q s q d (i)/l s (x), VxeR 2n 

S=1 

Qs h o, s = 0,... ,m. 

Pi G C, i = l,...,N 

Note that, as already discussed in Section 5.1, the constraint 
v 3 — coj x > 0 is satisfied for all x £ X. Therefore, the 
half-space: 7 '-Lj = {p€ R" : coj p < zz,-} is guaranteed to 
contain X. Thus, also the set M is included in 7 -Lj. Finally, 
note that, in order to deal with the nonconvex constraint 
LOj f 0 in (50), Problem (50) can be splitted into the two 


N 

LO* , v* = arg min R { n 0 } (ft) 

uje r 1=1 

Q s 

s.t. 

<-Oj 1 = —1 

T (51b) 

v :j - tOjX = q d (x) Q 0 q d (x)+ 

m 

- ^q d (x) T Q s q d (x)/i s (x), VxeR 2 ” 

S=1 

Q s h 0, s = 0,..., m. 

Pi € £, * = 1,... ,7V 

with iOj t i denoting the first component of vector to :] . The 
optimizer {u>*, v *} of Problem (50) is the given by the pair 
{uJ*,77*} or {to*,jg*} that provides the minimum value of 

the objective function EiLi R {Uj}{Vi)- 


Remark 5 For a fixed degree 2d of the SOS polynomials, 
the number of optimization variables of Problems (51) in¬ 
creases polynomially with the state dimension n and lin¬ 
early with the number m of constraints h s (x) defining the 
set X. Specifically, the number of optimization variables of 
Problem (51) is 0(mn 2d ). In fact, the number of free de¬ 
cision variables in the matrices Q s (with s = 0,..., m) is 


(2n+d 

\ d 


) (l + (“«)) 


= 0(n 2d ). On the other hand, for a 


fixed n, the size of the matrices Q s increases exponentially 
with the degree 2d of the SOS polynomials. In order not to 
obtain too conservative results, practical experience of the 
authors suggests to take d > [4] + 1, where [•] denotes 
the ceiling operator. We remind that d is the degree of the 
considered polynomial system in (1). Roughly speaking, be¬ 
cause of memory requirement issues, the relaxed SDP prob¬ 
lems (51) can be solved in commercial workstations and 
with general purpose SDP solvers like SeDuMi in case of 
polynomial systems with 4 state variables and of degree d 
not greater than 6. Systems with more state variables can 
be considered in case of smaller values of d. Similarly, sys¬ 
tems of higher degree can be considered in case of a smaller 
number of state variables. 
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Remark 6 A? already discussed, Algorithm 2 computes, at 
each iteration, an half-space Hj : u) ;l x — Vj <0 containing 
the set X (thus also A4), i.e., 

ujj x — Vj < 0 V x £ X. (52) 

The parameters U3j and Vj are then computed by solving 
Problem (50), and replacing the robust constraint (52) with 
a SOS constraint (see Problem (50)). Note that the same 
principles of Algorithm 2 and of the SOS-based relaxation 
discussed in this section can be used to compute, instead of 
an half-space TLj, a more complex semialgebraic set (e.g., 
an ellipsoid) described by the polynomial inequality: 

ut T q(it) <0 V x £ X, (53) 

with q(x.) being a vector of monomials in the variable x. The 
parameters u> can be then computed by properly modifying 
the SOS-relaxed Problem (50). For instance, in case we are 
interested in computing an ellipsoidal outer approximation 
of X, the function uj T q(x) should have a quadratic form, 
and its Tlessian should be enforced to be positive definite. 

Example 4 Let us continue with Example 2. Fig. 4 shows 
the polytope obtained by applying Algorithm 2 solving Prob¬ 
lems (51) instead of the nonconvex optimization in A2.2. 
The SDP Problems (51) are solved for a degree of the SOS 
polynomials equal to 2d = 4. The solution is a polytope 
S that outer-bounds X\. It can be observed that because 
of the approximations introduced (SOS and the approxima¬ 
tion of the indicator functions), which are necessary to ef¬ 
ficiently solve the optimizations, the half-spaces bounding 
X\ are not tangent to it and the computed region S still in¬ 
clude two black points. Therefore, the computed polytope is 
not the minimum-volume polytope. However, it is already a 
very good outer-approximation of it. In the next section, we 
describe a further refinement of Algorithm 2 aiming to com¬ 
puting a tighter polytope S. According to the steps Al.1.3 
and Al.1.4 of Algorithm 1, we outer-approximate AA (and 
so X\) with S. At the next time step (k = 2) of the set- 
membership filter, we repeat the procedure to compute a new 
polytope outer-bounding Xi- The difference is now that in¬ 
stead of hif) in (36), we have the 9 linear inequalities that 
define the polytope in Fig. 4. This procedure is repeated re¬ 
cursively in time. 

6.5 Refinement of the polytope S 

Summarizing, an approximate solution of the robust opti¬ 
mization problem (43) is computed by solving the convex 
SDP problems (51), and, on the basis of Algorithm 2, the 
polytopic-outer approximation S of the set AA is then de¬ 
fined as s = B n H 1 n... n Hj. 

Note that, in solving (51) instead of (43), two different 
sources of approximation are introduced: 

• Approximation of the indicator functions I i < y j i (pf) with 
the convex functions R^- H .y(pi) (see Fig. 3); 



Fig. 4. Final polytope after running Algorithm 2. 

• Approximation of the robust constraint v — uo x > 
0 Vx £ X with the convex conservative constraint 
v - u> T x = cr 0 (x) - cr s (x)/i s (x). 

The latter source of approximation can be reduced by in¬ 
creasing the degree 2d of the SOS polynomials. In fact, 
as already discussed in Section 5.1, according to the Puti- 
nar’s Positivstellensatz each function v — uj 1 x such that 
v — u; T x > 0 Vx £ X can be written as v — tu T x = 
rro(x) — cr s(x)/i s (x) provided that the degree of the 

SOS polynomials ag,ai,... ,a m is large enough. On the 
other hand, there is no theoretical result concerning the ac¬ 
curacy of the approximation of the indicator functions in 
Problem (43) with the convex functions R^ j y(pi) appear¬ 
ing in Problem (51). Because of that, the polytope S obtained 
by solving convex problems (51) (for j = 1 ,..., J) is not 
guaranteed to minimize the original nonconvex optimization 
problem (42). Algorithm 3 can then be used to refine the 
poly topic outer approximation S provided by Algorithm 2. 

The main principle of Algorithm 3 is to process, one by one, 
all the points belonging to the polytopic outer-approximation 
S initially given by Algorithm 2. For each of such points 
Pi, an half-space Hi : uj* x — a* <0 including the set 
X (i.e., X C Hi) and at the same not containing the point 
Pi (i.e., pi £ Hi, or equivalently — uj* p t + v* < 0) is 
seeked. In this way, all the points pi which do not belong to 
the minimum volume polytopic outer approximation of X 
are discarded. Thus, a tighter (but more complex) polytopic 
outer approximation of X is obtained. 

An important feature enjoyed by the refined polytope S* is 
given by the following theorem. 

Theorem 7 The polytope S* computed with Algorithm 3 is 
a global minimizer of problem (42). 

Proof: Let S be a polytope belonging to the set of feasibility 
of problem (42) (i.e., AA C S) which does not minimize 
(42). This means that there exists a polytope S such that 
M C S C S and a point p given as input of Algorithm 2 
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Algorithm 3: Refinement of the polytope S 
[input] Sequence of the random points pi provided as input 
of Algorithm 2 and such that pi £ S. Let N be the number 
of points pi belonging to S. 

A3.1 5* <- 5 
A3.2 for i = 1 : N 

A3.2.1 Compute the solution of the following optimiza¬ 
tion problem 

=arg min — uj T pt + v 

wer 


s.t. 

u) f 0 

v — u> T x. > 0 Vx £ X. 



■n 


Fig. 5. Exampe 1: hyperplanes defining the polytope <Sj (black 
lines) and true state uncertainty set X\ (gray region). 

problems (51) in the format required by the used SDP solver. 


A3.2.2 

[output] Polytope S*. 


such that: p £ S and p (f_ S. Thus, for p, = p, the optimal 
solution {uj* , v* } of Problem (54) is such that * Pi — V i > 
0. Let Hi be the half-space defined as Hi : u>* T x — v* < 0. 
Obviously, p Hi. Besides, the output S* of Algorithm 3 
is contained in the hyperspace Hi . Therefore, since p f 77., 
and S* C Hi, it follows that the point p S*. Then, a 
polytope S that does not minimize the optimization problem 
(42) can not be the output of Algorithm 3. ■ 

Theorem 7 mainly says that there exists no polytope includ¬ 
ing Ai and containing less randomly generated points p, 
than S*. However, it is worth remarking that only an ap¬ 
proximated solution of Problem (54) can be computed, as 
the robust constraint v — u)' x >0 Vx £ A appearing in 
(54) has to be handled with the SOS-based techniques de¬ 
scribed in the previous section. Thus, conservativeness could 
be added at this step. Therefore, the main interpretation to 
be given to Theorem 7 is that Algorithm 3 cancels the effect 
of approximating the indicator function /{- H }(p i ) with the 
convex function R^.y(pi). 

Example 5 Let us continue with Example 2. Fig. 5 shows 
the computed polytope Sf along with the true state uncer¬ 
tainty set X\. The CPU taken by the proposed algorithm to 
compute the 54 hyper-spaces that define the polytope 5* is 
about 830 seconds. However, only 80 out of 830 seconds 
are spent by the solver SeDuMi to solve 108 (i.e., 54 x 2) 
SDP problems of the type (51). The other 750 seconds are 
required by the SOStools interface to formulate, 108 times, 
the SDP problems (51 ) in the format used by SeDuMi. There¬ 
fore, the computational time required to compute the poly¬ 
tope SI can be drastically reduced not only by using more 
efficient SDP solvers, but also directly formulating the SDP 


7 Numerical examples 

Let us consider the discrete-time Lotka Volterra prey- 
predator model [63] described by the difference equations: 

x i( k ) = x i (k — 1) ( r +1 — r x i ( k — 1 ) — bx 2 ( k — 1 )) + w 1 (k — 1), 

X 2 (k)=cx\(k — l)a: 2 (fc — 1 ) + (1 — d)x 2 {k — l)+W 2 (k — 1 ), 

(55) 

where x±(k) and ^(fc) denote the prey and the predator 
population size, respectively. In the example, the following 
values of the parameters are considered: r = 0.25, b = 0.95, 
c = 1.1 and d = 0.55. The observed output is the sum of 
the population of the prey and predator densities, i.e., 

y(k) = x\(k) + x 2 ( k ) + v(fc), (56) 

where the measurement noise v(/c) is bounded and such that 
|| v(fc) ||oo < 0.05. The initial prey and predator sizes x(0) = 
[xi(0) tC 2 (0)] T are known to belong to the box Xq = 
[0.28 0.32] x [0.78 0.82] and the noise process w (k) = 
[wi(k) W2{k)] T is bounded by ||w(A;)Hoc < 0.001. The 
data are obtained by simulating the model with initial con¬ 
ditions £i(0) = 0.8 and £ 2 ( 0 ) = 0.3, and by corrupting 
the output observations with a random noise v(/c) uniformly 
distributed within the interval [—0.05 0.05]. 

Polytopic outer approximations Sf ; . of the state uncertainty 
sets Xk (with k = 1,..., 40) are computed through Algo¬ 
rithm 2. N = 20 random points are used to approximate the 
volume of the polytope S£ (as described in Section 6.1). In 
order to limit the complexity in the description of the poly¬ 
topes SI , the maximum number of halfspaces describing S* ; . 
is set to 8. This means that Algorithm 2 is stopped after at 
most 4 iterations (we remind that the initial outer-bounding 
box Bk is already described by 4 half-spaces). When the 
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Fig. 6. Example 2: outer-bounding polytopes (gray) and true state 
trajectory (black dots). 

output of Algorithm 2 is a polytope SJ. described by less 
than 8 half-spaces. Algorithms 3 is used to refine the poly¬ 
topic outer approximation Sf,. Fig. 6 shows the computed 
polytopes outer approximating the state uncertainty sets 
Xk (with k = 1, ..., 40), along with the true state trajectory. 
The Hybrid toolbox [64] has been used to plot the poly¬ 
topes in Fig. 6. The average CPU time required to compute 
a polytope S k is 28 seconds (not including the time required 
by the SOStools interface to formulate the SDP problems 
(51) in the format used by the solver SeDuMi ). For the sake 
of comparison. Fig. 7 shows the outer-bounding approxima¬ 
tions of the state uncertainty sets Xk when boxes, instead of 
polytopes, are propagated over time. For a better compari¬ 
son, in Fig. 8 the bounds on the time-trajectory of each state 
variable obtained by propagating boxes and polytopes are 
plotted. The obtained results show that, as expected, propa¬ 
gating polytopic uncertainty sets instead of boxes provides 
a more accurate state estimation. Finally, we would like to 
remark that a small uncertainty on the noise process is as¬ 
sumed (i.e., || w(fe) Hco < 0.001) since, for larger bounds on 
||w(fc)||oo, it would not be possible to clearly visualize the 
uncertainty boxes in Fig. 7. 

8 Conclusions 

In this paper we have shown that set-membership estima¬ 
tion can be equivalently formulated in a probabilistic set¬ 
ting by employing sets of probability measures. Inferences 
in set-membership estimation are thus carried out by com¬ 
puting expectations with respect to the updated set of proba¬ 
bility measures V, as in the probabilistic case, and they can 
be formulated as a semi-infinite linear programming prob¬ 
lem. We have further shown that, if the nonlinearities in the 
measurement and process equations are polynomial and if 
the bounding sets for initial state, process and measurement 
noises are described by polynomial inequalities, then an ap¬ 
proximation of this semi-infinite linear programming prob- 



Fig. 7. Example 2: outer-bounding boxes (gray) and true state 
trajectory (black dots). 



Fig. 8. Example 2: bounds on state trajectories obtained by prop¬ 
agating boxes (black line); bounds on state trajectories obtained 
by propagating polytopes (gray line); true state trajectory (black 
dots). 

lem can be obtained by using the theory of sum-of-squares 
polynomial optimization. We have finally derived a proce¬ 
dure to compute a polytopic outer-approximation of the true 
membership-set, by computing the minimum-volume poly¬ 
tope that outer-bounds the set that includes all the means 
computed with respect to V. It is worth remarking that the 
set-membership filtering approach discussed in the paper can 
be extended to handle noise-corrupted input signal obser¬ 
vations and uncertainty in the model parameters, provided 
that the corresponding state uncertainty set Xk remains a 
semi-algebraic set. As future works, we aim first to speed up 
the proposed state estimation algorithm in order to be able 
to use it in real-time applications in systems with fast dy¬ 
namics. To this aim, dedicated numerical algorithms, written 
in Fortran and C++, for solving the formulated SDP opti- 
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mization problems will be developed. Furthermore, the SDP 
problems will be directly formulated in the format required 
by the SDP solver, thus avoiding the use of interfaces like 
SOStools. An open source toolbox will be then released. 
Second, by exploiting the probabilistic interpretation of set- 
membership estimation, we plan to reformulate it using the 
theory of moments developed by Lasserre. This will allow 
us to ground totally set-membership estimation in the realm 
of the probabilistic setting, which will give us the possibility 
of combining the two approaches in order to obtain hybrid 
filters, i.e., filters that include both classical probabilistic un¬ 
certainties and set-membership uncertainties. 
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